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Abstract 

Due to their inherent dissipation and stability, the MacCormack scheme and its variants 
have been widely used in the computation of unsteady flow and acoustic problems. 
However, these schemes require many points per wavelength in order to propagate waves 
with a reasonable amount of accuracy. In this work, the linear wave propagation 
characteristics of MacCormack-type schemes are shown by solving several of the CAA 
Benchmark Problems. 


Introduction 

In the field of computational aeroacoustics, numerical schemes are expected to propagate 
waves accurately for long distances over long periods of time. In order to accomplish this 
goal, a certain number of spatial points are required per wavelength to model each wave, 
and a certain time step is required in order to model the wave's movement in time. It is 
desirable from a computational standpoint to reduce the number of points required per 
wavelength and increase the size of the allowable time step. 

One popular and well-tested method uses a modification of the MacCormack scheme 1 , 
which is second order accurate in time and fourth order accurate in space. This extension 
of the MacCormack scheme is known as the 2-4 scheme, and was described by Gottlieb 
and Turkel. 2 This scheme has been used successfully on a wide range of fluid and 



aeroacoustics problems. 315 Sankar, Reddy, and Hariharan have evaluated this scheme for 
aeroacoustics applications. 16 It has been extended to sixth- order spatial accuracy by 
Bayliss, et. al. (2-6 scheme) 17 , and an extension of the 2-4 scheme to fourth-order time 
accuracy is described by Viswanathan and Sankar. 18 

Building on this previous work, a new high-accuracy MacCormack-type scheme has been 
developed for use in computational aeroacoustics. 19 This scheme has been successfully 
applied to the real-world problem of supersonic jet noise prediction. 20 In this paper, the 
performance of this scheme will be evaluated using the benchmark problems of the first 
and second CAA workshops. The results are used to quantify the performance of the 
various schemes. 

Numerical Formulation 


In this work, four previously existing MacCormack-type schemes will be used: the 

classical MacCormack scheme, the 2-4 scheme of Gottlieb and Turkel, the 2-6 scheme of 
Bayliss, et. al., and the 4-4 scheme described by Viswanathan and Sankar. 

In addition to these schemes, the high-accuracy MacCormack-type scheme of Hixon will 
be used. This scheme utilizes the Dispersion Relation Preserving methodology of Tam 
and Webb 21 for the spatial discretization and the 4-6 Low-Dispersion and -Dissipation 
Runge-Kutta scheme of Hu, et. al. 22 for the time integration. It is formally fourth-order 
accurate in time and space for linear problems. 

The time marching method used by these MacCormack-type schemes can be written as 
follows: 

QO) = Q k 

Q (2) = Q k +a 2 A1F(q (1) ) 

Q< 3 >=Q* + a 3 AtF(Q< 2 >) 
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Q w =Q k +a 4 AfF(Q (3) ) 
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where the values of the coefficients are given in Table 1. Notice that the 4-6 time marching 
method alternately uses four and six stages to move to the next time level. Each spatial 
derivative uses biased differencing, either forward or backward, providing inherent 
dissipation for the solver. 

Using a spatial derivative at point j as an example. 


Forward: 


dQ 

dr 


A r 


a_iQ,_i + 3qQi 

+a A+1 + a 2®i+2 

+a 3^i+3 


( 2 ) 


Backward: 


dQ 

dr 


A r 


a 3@i-3 + a 2@i-2 

+a-|6/_i + a 0 Qj 

+a_iQ/ + i 


(3) 


The sweep directions are reversed between each stage of the time marching scheme to 
avoid biasing, and the first sweep direction in each time step is alternated as well. This 
gives a four-step time marching cycle (using the 4-6 time marching method to illustrate): 
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Table I: Coefficients for Runge-Kutta Time Stepping Schemes. 
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Table II: Coefficients for MacCormack-Type Schemes 
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Q k +' = L BFBF Q k 

Q k+2 =L FBFBFB Q k+ ' 

n* +3 - / n k+2 
U - l fbfb u 

n k+ 4 - / n*+ 3 

W - l bfbfbf u 

At the computational boundaries, flux quantities outside the boundaries are needed to 
compute the spatial derivatives; these fluxes are determined using third-order extrapolation 
from the interior nodes. 


Benchmark Problems 

These schemes are compared using linear benchmark problems from the first 23 and second 
CAA Workshops. 

1 ) One Dimensional Problems 

Problem 1 of the first workshop requests the solution at t = 400 of: 

u t + u x = 0 (5) 


where 


u(0) = -exp 


-20 < x < 450 



Ax = 1.0 


( 6 ) 


Problem 2 of the first workshop requests the solution at t = 400 of: 

u 

u t + u r + y = 0 


(7) 


where 
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u(r,0) = 0 
u(5,t) = sin( cot) 


n k n 

(o = 

3 4 6 
5 < r < 450 
Ar = 1.0 


( 8 ) 


These problems test the ability of the scheme to accurately propagate linear waves of 
varying wavelengths for long distances of travel. 

2) Two Dimensional Problems 


The problems given for the second CAA Workshop are much more difficult. In this 
problem set, a 2-D circular cylinder of radius 0.5 is placed at the origin, and acoustic waves 
reflect and scatter from this curved surface. The governing equations are the two- 
dimensional linearized Euler equations in polar coordinates: 
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For Problem 1, S is a simple harmonic source at r = 4, 0 = 0, given by: 


S = 


0 

0 


exp 


-ln(2) 


(x-4 f + y 2 ) ; 


( 0 - 2 ) 


sin(8xf) 


( 10 ) 


The problem requests rp 2 in the limit as r -> ®o in the arc 1 80 > 0 > 90. 


In Problem 2, S is an initial disturbance at time t = 0, given by: 
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For this problem, the pressure time history from 6 < t < 10 at three points is requested. 
These points are at 0 = 90° (point A), 0 = 135° (point B), and 0 = 180° (point C) at a 
radial distance of 5 from the origin. 

Boundary Conditions for Two-Dimensional Problems 


There are three boundary conditions which are used. At the cylinder surface (r = 0.5), the 
Thompson solid wall boundary condition is used, and the equations become: 
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( 12 ) 


In this computation, three ghost points are used inside the surface for the radial derivative; 
their values are set as: 
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In the far field (r = R max ), the acoustic radiation condition is used: 
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(14) 


For the radial derivative at the outer boundary, three ghost points are used. The values of 
the variables at these ghost points are determined using third-order extrapolation from the 
interior values. 
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At the symmetry planes (0 = 0 and 0 = 7t), a symmetry condition is used. For example, 
around i = 1 : 
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V 1 


> = - 

-Vb‘ 

.P‘, 

1-1 

. ft . 


where i is the index in the azimuthal direction. 


(15) 


Computational Grid for Two-Dimensional Problems 

For Problem 1, a 801 (radial) x 501 (azimuthal) grid was used, covering a domain of 0.5 < 
r < 20.5 in the radial direction, and 0 < 0 < tt. Since the wavelength of the disturbance is 
0.25, this grid results in 7-10 points per wavelength. The exact results were given at the r = 
15 line, giving a maximum of 76 wavelengths of travel at 0 = n. 

For Problem 2, a 201 (radial) x 301 (azimuthal) grid was used, covering a domain of 0.5 < 
r< 10.5 in the radial direction, and 0 < 0 < n. Since the transient problem only requires 
data from 6 < t < 10, the outer radial boundary only has to be far enough away such that no 
reflections can reach any of the three data points during this time period. 

Results 


11 One Dimensional Results 

The solutions of Benchmark Problem 1 at t = 400 as calculated by the various schemes are 
shown in Figures 1-5. In this problem, a Gaussian pulse propagates in time and space for 
a given time, at which point the results are compared. On each figure, the solution is 
shown for various time steps to illustrate the effect of the time step on the accuracy of the 
time integration. 

Figure 1 shows the solution obtained using the classical MacCormack scheme. Two 
points are evident: first, the classical scheme is a perfect propagator at a CFL number of 
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1.0; second, the scheme is very dissipative and dispersive for other time steps. Since 
perfectly uniform grids are unusual in real-world problems, the perfect propagation is not 
overly useful. 

Figure 2 shows the solution obtained using the 2-4 scheme of Gottlieb and Turkel with 
time steps of 0.25, 0.4, and 0.5. The dissipative nature of the scheme is evident, and the 
time step has a large effect on the dispersion error. 

Figure 3 shows the solution obtained by the 2-6 method of Bayliss, et. al. Again, the 
solution shows dissipation and dispersion, with higher dispersion errors than the 2-4 
scheme. 

Figure 4 shows the solution obtained by the 4-4 method described by Viswanathan and 
Sankar. The effect of the increased time accuracy on the dispersion error is immediately 
apparent. Another point of interest is that the four-stage, fourth-order Runge-Kutta scheme 
is dissipative at larger time steps. 

Figure 5 shows the same time-stepping scheme using the optimized DRP spatial 
differencing. The solution obtained is very accurate and is also insensitive to the time step 
chosen. 

Figures 6-8 compare the results of the 2-4 scheme to those of the optimized DRP scheme 
for Benchmark Problem 2. In this problem, a single-frequency spherical wave propagates 
outward from an impulsively-started vibrating sphere. As given, the problem has two 
parts: the first has 8 points per wavelength, and the second has six points per wavelength. 
To give the 2-4 scheme some chance of obtaining reasonable results, a third part was added 
which has 12 points per wavelength. 

Figure 6 shows the results for 12 points per wavelength at 33.33 wavelengths of travel. 
The solution given by 2-4 scheme shows dissipation error and some dispersion error, 
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while the solution given by the optimized DRP scheme shows very little dispersion or 
dissipation error, even with a time step five times larger than that taken by the 2-4 scheme. 

Figure 7 shows the results for 8 points per wavelength at 50 wavelengths of travel. The 
solution given by the 2-4 scheme is very dissipated and shows some dispersion error. The 
solution given by the optimized DRP scheme shows some dispersion and dissipation 
error, but is still very good. 

Figure 8 shows the results for 6 points per wavelength at 66.67 wavelengths of travel. 
Except for the initial transient, the wave has been completely damped by the 2-4 scheme, 
while the optimized DRP scheme is still showing good accuracy in dispersion. However 
the wave has dissipated to approximately 65% of the exact value. 

2) Two Dimensional Results 

Results for Problem 1 of the second CAA Workshop are given in Figure 9. In order to 
avoid problems with the very large initial transient, a polynomial function was used to 
smoothly increase the amplitude of the forcing function. The time step used was limited 
by the stability of the solid wall boundary; for these calculations a time step of At = 
0.00245 was used (CFL = 0.786). The calculation was run to a time of 32.09, with results 
being taken from 31.59 < t < 32.09. This calculation took a total of 6.27 hours of CPU 
time on a Cray Y/MP, running at 191 Mflops. The results are given at r = 15 D, and 
compare very well with the exact solution. 

Results for Problem 2 are given in Figures 10-12. The results agree very well with the 
exact solution. This calculation, using a At of 0.0025 in order to print out the required 
results, took a total of 469 CPU seconds on a Cray Y/MP, running at 175.5 Mflops. 
However, the code could run stably at a CFL number of 0.864 (At = 0.0045), requiring 
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261 CPU seconds. With more stable solid wall boundary conditions, it is expected that the 
scheme can recover the CFL = 1 .4 time step that has been seen previously. 

Grid refinement studies were conducted for Problem 2; the effect of halving and doubling 
the grid are shown for Point C in Figures 13 and 14. Point C was chosen because it was 
the most distant point from the initial location of the pulse. In Figure 13, three computed 
results are shown: a half grid (101 x 151), the grid used (201 x 301), and a doubled grid 
(401 x 601). The two denser grids have nearly identical results, and compare very well 
with the exact solution. The coarsest grid, however, shows leading and trailing waves, 
some traveling much faster than the physical wave. This is due to the low resolution of the 
grid causing the solver to incorrectly allow high-frequency waves to travel faster than the 
speed of sound. 

Figure 14 shows the transient peak at point C. The effect of increased grid is illustrated in 
this graph; the transient peak becomes closer and closer to the exact solution as the grid 
becomes denser. At this extreme amplification, it can be seen that the transient peak 
velocity is very slightly off with the grid used, but the answer is well within expected 
tolerances for this case. 

Conclusions 

In this work, the dispersive and dissipative characteristics of a new high-accuracy 
MacCormack-type scheme were investigated using benchmark problems of the first and 
second CAA workshops. The results show that this new scheme is very promising for 
computational aeroacoustics applications, requiring only 6-8 points per wavelength and 
allowing large time steps. 
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The MacCormack-type schemes are of great interest due to their ease of programming and 
use, and inherent numerical dissipation. This work shows that this type of scheme can be 
optimized to perform very well. 
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Figure 1 . — Solution of Benchmark Problem 1 of the First CAA 
Workshop using Classical MacCormack Scheme. 
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Figure 2. — Solution of Benchmark Problem 1 of the First CAA 
Workshop using 2-4 Scheme. 
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Figure 3. — Solution of Benchmark Problem 1 of the First CAA 
Workshop using 2-6 Scheme. 
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Figure 4. — Solution of Benchmark Problem 1 of the First CAA 
Workshop using 4-4 Scheme. 
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Figure 5. — Solution of Benchmark Problem 1 of the First CAA 
Workshop using Optimized DRP Scheme. 
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Figure 6. — Comparison of Solutions of Benchmark Problem 2 
(12 points/wavelength). 
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Figure 7. — Comparison of Solutions of Benchmark Problem 2 
(8 points/wavelength). 
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Figure 8. — Comparison of Solutions of Benchmark Problem 2 
(6 points/wavelength). 
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Figure 9. — Computed solution of Problem 1 of the Second CAA 
Workshop at r = 15 using Optimized DRP Scheme. 
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Figure 10. — Computed solution of Problem 2 of the Second CAA 
Workshop at point A using Optimized DRP Scheme. 



Figure 1 1 . — Computed solution of Problem 2 of the Second CAA 
Workshop at point B using Optimized DRP Scheme. 
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Figure 12. — Computed solution of Problem 2 of the Second 
CAA Workshop at point C using Optimized DRP Scheme. 



Figure 13. — Effect of grid density on computed solution at 
point C. 
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Figure 14. — Effect of grid density on computed solution for 
transient peak at point C. 
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